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Abstract 

Powell’s nonlinear programming code, VF02AD, has been used to generate approximate 
minimum-time tip trajectories for 2-link semi-rigid and flexible manipulator movements in the 
horizontal plane. The manipulator is modeled with an efficient finite-element scheme for an n-link, 
m-joint system with horizontal- plane bending only. Constraints on the trajectory include bound- 
ary conditions on position and energy for a rest-to-rest maneuver, straight-line tracking between 
boundary positions, and motor torque limits. Trajectory comparisons utilize a change in the link 
stiffness, El , to transition from the semi-rigid to flexible case. Results show the level of compliance 
necessary to excite significant modal behavior. Quiescence of the final configuration is examined 
with the finite-element model. 

Introduction 

Trajectory planning is essential in budgeting a manipulator’s actuator efforts to maximize 
productivity. For repetitive tasks, the minimum- time maneuver goes hand-in-hand with this goal. 
A variety of approaches have been advanced for rigid manipulator control, taking advantage of the 
fact that all or some of the controls take the form of switching functions between actuator bounds. 
Bobrow [1] used an intuitive approach to generate optimal switching controls, as well as proving 
the boundedness of the controls. Weinreb and Meier [2] [3] used calculus of variations approaches 
to incorporate control bounds in the problem formulation. In a second study, Bobrow [4] used 
numerical optimization to generate spline fits to the switching controls. 

Switching functions do not lend themselves to maintaining tip accuracy for non-rigid struc- 
tures. One would hope that the applied controls do take advantage of the bounds to maximize 
performance, but a clear analytical directive for this does not exist at the present. 

In filling this void, parameter optimization techniques can provide approximate optimal perfor- 
mance solutions for systems driven by complex, highly nonlinear dynamic models with arbitrary 
equality or inequality constraints. Of these solution techniques, Powell’s Recursive Quadratic Pro- 
gramming algorithm [5], embodied in the code VF02AD, has proven to be a robust tool for a variety 
of aerospace applications [6] [7] [8], and will be used in this study. The primary drawback to this or 
other numerical optimization methods is the dependancy on accurate gradient approximations of 
the performance index and constraints with respect to the parameters. 

The ensuing discussion initially describes the structural dynamics model of the manipulator, 
followed by the optimal control problem and parameterization of the controls, and ends with the 
results of a computational experiment. 

The manipulator structure modeled in this study, and presently under fabrication, is a 2-link 
cantilever arrangement constrained to slew in the horizontal plane. Tall, thin links are used to 
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minimize vertical plane droop. The hub or joint- 1 actuator slews both links, an interlink motor, 
and tip payload. The interlink or joint-2 actuator located at the end of link-1 slews the second 
link and the tip payload. The joint- 1 /joint-2 actuator torque ratio is about 4/1. The complete 
manipulator is about 0.5 meters (m) tall and 1.2m long. 

The Structural Model 

There is a long literature discussing the difficulties of simulating the vibrations of rotating 
structuxes[9] [10] [11]. The problem seems to arise from kinematics that are of second order impor- 
tance in nonrotating problems, but become of first order importance in the presence of rotational 
accelerations. Additonally, there are constraints inherent to the flexible link problem which must 
be satisfied: motions occur entirely in a horizontal plane; one end of the chain of links is attached 
to a stationary hub; and each flexible link is inextensible. These considerations motivate the de- 
velopment of a mathematical model that faithfully carries the full kinematics of the problem. It is 
also necessary to devise such a model in a form that will lend itself to real-time calculations. 

The need to meet these apparently conflicting demands motivated the development of a model 
specialized to flexible, multilink structures. That apparently successful strategy is outlined below. 
The full kinematics are retained by expressing the configuration as functions of convected coordi- 
nates. This is a traditional approach in nonlinear elasticitity [12]. Further, the kinematic variables 
are selected so that all geometric constraints (fixed hub, planar motion, and non-extension) are 
automatically satisfied. 

Since motions are assumed to occur entirely in a plane, it is also assumed that the elastic lines 
of the links as well as the mass centers of the cross sections all lie in the same plane. Each cross 
section is identified by its arc-length distance from the hub, so that the orientation of the center of 
the cross section s at time t is 


fi(s f t) = cos(0(3,t))x + sin(0(s,f))J 


The location of the center of cross section $ at time t is obtained by integration of the above unit 
tangent vector: 


*(M) = / t)da' 

Jo 


Similarly, the velocity at the cross section s at time t is obtained by integration of the time derivative 
of 0(a, t): 


x(a,t) = f 6(a\ t)y(s\ t)da' 
Jo 


where 

7 (s, t) = — sin(0(s, t))i + cos(0(s, t))j 


The above description of configuration - entirely in terms of 0(s, t) - causes all of the geometrical 
constraints to be satisfied automatically. Additionally, the above description expresses the configu- 
ration in terms of one unknown field (0), instead of the more conventional two or three fields (®, t/, 
& 0 ). 

The governing equations of the dynamics are derived using those kinematics and a frame- 
invariant variational method - Hamilton’s principle. A finite element discretization is used to cast 
the resulting integro- differential equations for 0(s, t) and its first and second derivatives into a sys- 
tem of fully- coupled, nonlinear algebraic equations. Particularly important for the application at 
hand is the observation that since all spatial integrals are with respect to the convected coordinate, 
s, those integrals are configuration-independent and need be done only once. The nonlinearities 
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remain, and a new nonlinear system must be solved at each time step, but the time consum- 
ing quadrature process can be done in advance of the dynamics simulation. These features are 
illuminated below through derivation. 

Hamilton’s principle is that 

6 f h [j KE(t ') - SE(t') + WE{t')\ dt ' = 0 (1) 

Jti 

where KE is kinetic energy, SE is strain energy, and WE is external work. Quantities associated 
with the end times, and t 2 have to do with initial and final conditions and the conservation 
of momenta, but the governing equations necessary to model motion of the flexible structure are 
obtained by consideration of the integrand alone: 

SKE(t) - 6SE{t) + 6WE(t) = 0 (2) 


for all t\ < t < < 2 - 

The kinetic energy is that of the flexible links plus that of all concentrated masses and concen- 
trated moments of inertia: 


1 m ij < mua &cs . inertias 

KE i t )=~f P( a )*(M) -x(a,t)da + - ^ M k z(a k ,t) • x(a k ,t) + - l t 0{ai,t) • 0(a h t) 

ZJo 1 k=i 1 l=i 


k=i 

The strain energy is that of the flexible links 

rL 90(a,t) d0(a,t) 


sE{t)= \L n{a) 


ds 


ds 


ds 


Discretization of the above energy terms is obtained by discretizing the tangent vector (3 as: 

nodes 

Pi 3 ’ 1 ) = E Pri{t)Pn{a) 

n= 1 

where the shape functions, p n , have support over intervals that are small relative to the anticipated 
radii of curvature. The above condition on the support of the basis functions is necessary to assure 
compliance with the condition of nonextension. The resulting energies are: 


1 nodes nodes 

KE{t) = - Y, E 9m{t)0n{t)1mit)-Ut)Mm,« 


m=l n=l 


and 


where 


nodes nodes 


®(«)=» E E Pmit ) ■ fat) K m , n 


( 3 ) 


( 4 ) 


m=l n=l 


m ij masses meruaj 

M m,n = / p{s) q m (s) q n (s)ds + E Mk 9m{ 3 k) q n {sk) + E I l^K{l,rn)6 K {l,n) , 
Jo k = l l=i 

q m {a) = f p m (a)da , K mt7l = f n(a) p' m (a) p! n (a)da , 

Jo Jo 

6jc is the Kronecker delta function, i is a dummy variable, and p f m n are spatial derivatives. 


mertta« 
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After appropriate integration by parts, the integrand of equation 2 becomes: 

Tlodc9 

E Mm(— ?m(t) ■ P n (t)M m , n - 7m(0 • Pn(t)K m ,n + 7m(0 * *»( ) = 0 (5) 

n=l 

for all nodes m. In the above equation, r n is the external torque applied at node n. After /3 n (t) is 
expanded: 

= k(tmt) + («n('))V„(<) 

and Equation 5 is invoked for all a complete set of nodes second order equations in the nodes 
unknowns, 0 n , results as follows: 

nodes 

E 7m(0 • in(t)e n (t)M m ,n ~ f, m(<) * + 7m(*) • &(0*m.» - fm(«) * *m(«) (6) 

n=l 

The above problem formulation, involving only one unknown field, automatically satisfying all 
constraints, and requiring only one evaluation of element mass and stiffness matrices, lends itself 
to rapid numerical calculation. A computer code to generate and solve the above described system 
of equations for each time step has been written, tested and used by VF02AD. Both the derivation 
and code mentioned above are described more fully in Ref. [13]. 

One particularly interesting calculation, discussed in the literature as being difficult to solve, is 
that of a beam accelerated around a hub to an angular frequency above the first bending frequency 
of the beam. This is a particularly stringent test of flexible-dynamics codes, testing numerical 
robustness, as well as the correctness of the physics. Figure 1 shows the motion of the tip of such a 
flexible beam relative to the tip of a rigid beam rotating at the hub velocity. Initially, the flexible 
beam lags its rigid counterpart; it then overtakes and oscillates about the rigid beam. These results 
are in near exact agreement with the results presented in [9], obtained from a much more complex 
model. 

Optimal Trajectory Shaping 

The principle goal in this study is to combine the physics of the structure with optimiza- 
tion techniques to generate actuator torque histories for accomplishing a useful task with minimal 
degradation in performance. A secondary objective is to shortcut the work of an eventual feedback 
controller, which will be needed to compensate for modeling errors. 

Looking towards maximizing productivity in some repetitive task, a minimum-time tip tra- 
jectory was chosen for investigation. Constraints on such a trajectory include: completing a 
rest-to-rest maneuver, tracking a specified path (jr(f),y(t))et JM slewing between specified endpoints 
[(*(*<>)> y(*o)), (*(t/),y(t/))Jt» pi and not exceeding actuator torque limits r lt2niam - 

The configuration initially starts at rest. Driving a flexible structure to rest at the final time, 
tf y necessitates end constraints on both kinetic and potential or strain energies (KE(tf),SE(t/)). 
The chosen path is a straight line and actuator torques limits are constants. The problem can be 
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restated as 

minimize: J = tf 

- finite element model 

subject to: - input actuator torques, Ti |2 (t) 

- known initial conditions 
constrained by: 


■ 


Ztip(tf) - Z»pec«/«ed(*/) 

— 

' ■ 



Stap(f/) - y,pedfud{tj) 





Io f [ytip(x t ip(t)) - yw(* t <p(0)] rft 

= 

0 


= 

KE{t f ) 

= 




SE(t/) 

= 




Jo'(M*)l - 

< 


. 

j= 1,7 


< 

_ 


Note that the equality tracking constraint, C3, and inequality torque constraints, Ce,C7, are for- 
mulated as integrals. In addition, equality constraints on energy are point constraints. Both of 
these items will have profound effects on the example trajectories to be generated. 

Parameterization of the Controls 

To approximate optimum system performance from the aforementioned structural model, a 
suitable parameterization of the controls, Ti ( 2 (£), is necessary. The choice of “parameterizable ” 
torque functions is essentially limitless. For this study, the simplest case of using tabular values of 
torque, r t , as parameters at equal-spaced fixed times, for both joints was chosen, or 

n (ti) 7 T 2 (ti), * = l,n 0 <ti<tf, 

which results in 2 n control parameters. 

However, since the final time is changing due to minimization, the loss of control history def- 
inition would result if the times at which the control parameters are defined remain fixed in an 
absolute sense. To correct this, T it2 were specified at equally-spaced, nondimensional node points, 
C* = ti/tfj where 

n(Ct)iT 2 (C.), i=l,n 0<Ci<l, 

This allows the torque histories to “stretch” naturally over the trajectory length. Using this mod- 
ification, it is necessary to add tf as a parameter also, resulting in 2n + 1 control parameters to be 
found. 

Numerical derivatives of the performance index, tf , and the constraints, Cj(tf ), provided to 
VF02AD are central finite-difference approximations. To generate derivatives with respect to a 
torque parameter, r\ %2i = ti j2 (^), the parameter is perturbed equidistant about its nominal value, 
and complete trajectories (or integrations of the structural equations) are computed to the current 
nominal tf to produce perturbed Cj(tf) values. Since derivatives are computed over the current 
fixed tf , the derivatives, d£//dr 1>2 - = 0, and only the derivatives, 9Cj(tf)f dri f2i ^ 0. Naturally 
both tfjCj(tf) gradients with respect to £/, evaluated over the current nominal torque histories, 
are nonzero, (where dtf/dtf ~ 1.). 

Results 

The following finite-element structural model was used for the manipulator to produce the 
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sample trajectories. 


ITEM 

COMPOSITION 

LENGTH 

(m) 

MASS 

( k g) 

STIFFNESS, El 
(newton-m 2 ) 

joint ! bracket 

1 element 

.0635 

.545 

" 10 s 

link-1 

3 elements 

.5040 

.640 

10 2 ,10 3 

1st joint-2 bracket + 

1 element + 

.1070 

5.415 

10 6 

joint- 2 

2nd joint-2 bracket 

1 point mass 
1 element 

.1040 

.830 

10 5 

link 2 

3 elements 

.4890 

.313 

10 2 ,10 3 

Totals: 

9 elements 

1.2675 

7.743 



The two values of stiffness, j EL, for links 1,2 represent the trajectory comparison for this study. 
The brackets, modeled with a stiffness of 10 5 , are considered rigid. Point moments of inertia were 
used to define mass distribution for the brackets. No payload was used in this comparison. The 
joints were assumed to have no compliance or damping. 

The two trajectories, computed on a CRAY-XMP, were integrated for 100 time steps, where 
At = .01 tf. Trajectory evaluations for gradient computations executed in 0.75 secs. The torque 
histories for each joint were composed of 21 tabular values, where — .05. Torque bounds were 
chosen as ±16, ±4 newton-m (n-m) for joints 1 and 2. The path to be tracked for this study was 
the line connecting (x,y) pairs, (0.0,1.13) and (1.13,0.0). A composite of the slew motion for the 
“flexible” case ( EIi inkt = 100 newtons-m 2 ) is given in Fig.2. The parameterized torque histories 
that created this slew represent 100 iterations of VF02AD after initialization with the parameter 
solution values from the “semi-rigid” case ( EIunk» = 1000 n-m 2 ). The tip path traced is reasonably 
straight, but does contain some small ripples - a result of the integral statement of the tracking 
constraint, 

Figure 3 graphically depicts the difference between the semi-rigid and flexible links. Shown 
is the angular velocity of the finite-element node adjacent to joint 1. The frequency of vibration 
for Elunkg = 100 is about 19 hz. From examination of Fast Fourier transforms (FFTs) of the 
finite-element output of the system disturbed about the initial, midtime, and final positions, this 
appears to be one of the lower modes. Note the low angular velocity of the semi-rigid system after 
t/tf = 0.9, implying that a significant amount of time is being expended in order to bring the 
system to “rest”. This phenomena is definitely at odds with purely rigid system behavior. The 
KE(tf) = 0 constraint imposes a zero final angular velocity in both cases. Also, note that tf for the 
flexible case is slightly greater than the semi-rigid one. This demonstrates the approximate nature 
of the solutions, insofar as the semi-rigid solution not converging as well. 

Fig.4 shows the tj profiles. These still retain some of the boundedness qualities of purely rigid 
configurations. However, they begin and end near zero instead of the bounds (±16 n-m), and the 
transition between bounds is comparitively gradual. The oscillatory behavior in the Elunkg = 190 
torque is probably counteracting the excitement of the lowest structural modes, which would have 
the greatest Impact on the position constraints. Note the abruptness of the controls near the end 
in an attempt to quiet the structure. The r 3 torques in Fig.5 show minimal activity for most of 
the trajectory, except close to the end in order to accomplish the rest state. Note that this closing 
maneuver starts sooner with the more flexible link structure. 

The straight-line tracking error in millimeters (mm) is shown in Fig.6. Both torque histories 
appear to limit the error to ±5 mm except near the end of the EIu n ks — 1900 trajectory, where 
the error momentarily “escapes” before returning to zero to satisfy the end position constraints, 
(Ci(«/), C 3 (t/)). One drawback to the integral formulation is that it can relax tracking performance 
in isolated parts of the trajectory, yet yield a reasonably low residual (« 0) for C${tf). It may be 
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necessary to add interior point constraints to significantly decrease this error. A less complex 
alternative is to reduce the EIu n k 3 more gradually between solutions, in the fashion of a homotopy 
scheme, until the desired value is reached. 

The kinetic energy of the flexible structure at peak rates (t/tf « 0.45) is, not surprisingly, 
higher than the semi-rigid one as shown in Fig. 7, It is interesting that KE appears to be devoid of 
oscillatory behavior in both cases. Note again, the rest phase of the trajectories above t/tf = 0.9. 
Strain energy is shown in Fig. 8. This again displays the contrast first seen in Fig. 3. The semi-rigid 
structure produces relatively little strain, while the flexible-link configuration again contains the 
19 hz mode with a sizable increase in energy magnitude. Note the peaks in both cases, mirroring 
the sharp r r changes in Fig. 4. Also, note the enforcement of the SE(tf) = 0 point constraint at the 
end. The SE “floor” in both cases corresponds to the KE peaks in Fig. 7. 

The finite-element model was also able to examine the quiescence of the final state. At joint 
node accelerations can be zeroed and joint node positions fixed at their nominal final values. The 
corresponding torques applied to maintain these values can then be computed. Fig. 9 shows the 
residual T\ being applied after the approximate optimal tf = 1.505 secs for EIh^b = 100 has been 
reached. The oscillations in this residual appear to contain frequencies in the vicinity of 2.5 and 35 
hz. The lower appears to correspond to a fixed-free mode of the first link with the remainder of the 
system mass concentrated at the end. This mode was not seen during the optimized part of this 
trajectory. The second frequency may be a higher harmonic of the mode seen earlier, or possibly 
a numerical artifact. The small magnitude of this residual coupled with the spacing of the modes 
may be ideal targets for attack by a linear feedback control scheme. As an open-loop alternative, 
augmenting the tabular torque controls with frequency-based, parameterized functionals may be 
sufficient to suppress these oscillations. 

Conclusions 

A robust, parameter optimization tool has been able to generate actuator torque histories 
for approximate, minimum- time slewing maneuvers containing a variety of continuous and point 
constraints for a 2-link flexible manipulator. The parameters, or actuator torques, for each link 
were tabular values at fixed node points during the maneuver. Perturbations were made to each 
parameter to approximate final time and constraint gradients. The efficient formulation of the 
finite-element model made the numerical optimization procedure a realistic endeavor. 

The accuracy of the straight-line tip tracking was good. Additional interior constraints and/or 
a more gradual change in the stiffness should yield further improvements. For the trajectory used 
in this study, joint-1 applied most of the input, which included cancelling lower-mode vibrations 
for the structural as a whole. Energy constraints were effective in bringing the structure to rest at 
tf. It was also demonstrated that final energy constraints do not preclude vibrations during the 
slew. The intended production use of the manipulator will dictate whether this is a hindrance or 
not. Finally, the secondary goal of providing enough vibration control to maximize the success of 
a linear feedback controller in treating residual oscillations appears feasible. 
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Figure 3: Joint 1 node angular velocity 


D=1000. 

TF = 1514 SEC 

0=100. 

TF = 1505 SEC 



NORMALIZED TIME 1/1 F 

Figure 4: Tj joint torque histories 























